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FLUTTER ANALYSIS OF FLAT RECTANGULAR PANELS 
BASED ON THREE-DIMENSIONAL SUPERSONIC 
UNSTEADY POTENTIAL FLOW 

By Herbert J. Cunningham 
Langley Research Center 

SUMMARY 

A systematic analytical procedure has been developed for computing flutter char- 
acteristics of rectangular panels with stream-alined side edges, based on air forces 
from three-dimensional linearized supersonic unsteady potential flow. The procedure 
has particular usefulness in the low supersonic speed range where static and quasi- 
static aerodynamic approximations are considered to be least valid and can provide 
bases of comparison for some of the simpler types of analysis. The panel is considered 
to be finely divided into many boxes, and the aerodynamic influence coefficients between 
all pairs of boxes are obtained by numerical integration. The flutter analysis is a modal 
type, which readily coordinates with the aerodynamic box method, and can be used for 
calculating the flutter stability of any flat or nearly flat panel, whether of isotropic or 
anisotropic stiffness, and of buckled panels for which the flutter is a small- amplitude, 
simple harmonic, superimposed motion to which linear theory is applicable. A number 
of results are presented for flat unstressed, isotropic panels with simply supported edges 
and with clamped edges. For clamped-edge aluminum panels with a length-width ratio of 
2 at sea level, the panel flutter parameters are tabulated for eight Mach numbers ranging 
from 1.02 to 2.0. For Mach 1.3, flutter boundaries are plotted for length-width ratios 
from 0 to 10 for simply supported edges and from 0 to 4 for clamped edges so that 
design values can be read for a wide range of panel materials and air densities. Appen- 
dix A delineates the way in which the natural mode characteristics were developed for 
calculating the presented flutter results without the need for double-precision arithmetic. 
Appendix B provides formulas for conversion among a number of types of flutter solution 
parameters in current use. Appendix C describes a way to economize computer time 
for the large matrix multiplications required. 

INTRODUCTION 

A description of an analysis of the flutter of flat rectangular unyawed panels on the 
basis of unsteady potential flow and some results from that analysis were presented in 
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reference 1. Since its publication, considerably more information and results have been 
obtained. The purpose of this report is to present the overall results. The description 
of the analysis is sufficiently complete to enable the analyst to understand it without 
referring to reference 1. 

In the lower supersonic Mach number range, below a Mach number of about ^2, 
simple approximations for the aerodynamic forces, such as Ackeret theory or piston 
theory, do not give satisfactory results, at least for panel length-width ratios less than 
about 1.0. For length-width ratios greater than 1.0, there is now some evidence that 
Ackeret and piston theory can give good estimates of the panel thickness required to pre- 
vent flutter, and this evidence is discussed in the present report. 

The central feature of this analysis is the way in which the aerodynamic forces are 
determined. In what has come to be known in the literature as a "box method," the 
required double surface integration is carried out by considering the panel to be finely 
divided into a large number of equal- size rectangular elements or "boxes." The aerody- 
namic influence coefficient relating the velocity potential of each box to the motion of 
each other box is calculated by numerical integration and used in a modal-type flutter 
stability solution. Mode- shape properties from either calculation or experimental meas- 
urement can be applied in the flutter solution. 

A variety of flutter results are presented and discussed for flat unstressed panels 
with isotropic stiffness. Calculated mode shapes and frequencies were used. Although 
the results presented are limited to simple panels, this type of modal analysis with a box 
method is broadly applicable to flat or nearly flat panels, whether unstressed or stressed 
(as by thermal expansion), and whether of isotropic or anisotropic stiffness. It would 
also seem to be applicable to buckled panels for which the flutter is a small- amplitude, 
simple harmonic, superimposed motion to which linear theory would apply. 

A brief discussion is given of recently published findings that the stiffness coupling 
effects between assumed beam natural modes, when they are applied to the flutter analysis 
of clamped- edge plates, are not uniformly negligible and can become significantly large 
as length-width ratio increases. 

Appendix A presents the expressions for the panel vibration mode shapes and fre- 
quencies used in the flutter analysis in such forms that only single-precision arithmetic 
is necessary for all modes used. Appendix B gives the relations among some panel 
flutter parameters that are used in the literature. Appendix C describes a way to econ- 
omize computer time for the large matrix multiplication in computing the velocity 
potentials. 
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SYMBOLS 

a s speed of sound 

Ay aerodynamic influence coefficient giving the velocity potential at a box induced 

due to unit downwash on another box (eqs. (13) and (18)) 

Bs> B xs number of boxes in stream and cross- stream directions, respectively 

D flexural stiffness of an isotropic panel 

E Young's modulus of elasticity 

g,gj coefficients of mechanical hysteretic structural damping 

gj = Si + g 

hj,hj amplitude of natural mode- shape deflection for modes i and j 

Hj(x,y, r) distribution of time-varying panel deflection for mode j (eq. (2)), positive 
downward 

I G1 ,I G2’ I G3 integrals in Ay (eqs. (19) to (21)) 
k e reduced frequency with reference length e , ^ 

k 7 reduced frequency with reference length l , ^ 

6 y 

Kp>Kq roots of characteristic equations (eqs. (A3), (A7), and (A13)) 

^Pjq nondimensional eigenvalue quantity for natural vibration modes (eqs. (A4), 

(A12), (A17), and (A18)) 

l length of panel in stream direction 

m^ mass of panel per unit of surface area 
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M 


Mi 


Mach number, 

’ a s 

generalized mass for mode i 


* 

Mi = 7777^7 (applicable for uniform panel only) 


Zw m A 


p,q 


number of half-waves in panel vibration mode in stream and cross-stream 
directions, respectively 


Ap 


perturbation lifting pressure, positive downward 


Api 


q 


perturbation pressure due to mode j (eqs. (5) and (7)), positive sense same 
as that for Hj 

0 y2 

free- stream dynamic pressure, — 


qi(r),qj(r) time-varying generalized coordinate of motion for modes i and j, 
respectively 




complex amplitude of q^r) (eq. ( 2 )); also the eigenvector for mode j from 
a flutter solution 


Q 


i] 


generalized aerodynamic force from pressure due to mode j and modal 
deflection of mode i (eqs. ( 8 ) to ( 10 )) 


% 

S 


nondimensional computational quantity contained in Qyj (eq. (28)) 


panel area 


t thickness of panel 

u,v transformed panel coordinates and variables of integration in x- and 

y-directions, respectively, based on e/2 as a reference length, (eq. (15)) 

Uj,U 2 ,Vj,V 2 lower and upper limits of integration with respect to u and to v 
(eq. (16) and table I) 


speed of undisturbed airstream 
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w 


w i 


x,y 


x e ’?e 


x m>yn 


x,y 


X jjYj 


z 


a 



a* = 



a p,x’°q,y 


width of panel in cross-stream direction 

downwash velocity at panel surface for mode j, positive downward 

panel coordinates in stream and cross-stream directions, respectively 
(see fig. 1) 

panel coordinates based on e as a reference length, x e = *, Y e 

values of x e and y e at center of box m,n (see fig. 1) 

panel coordinates based on reference lengths l and w, respectively 

x- and y-variation, respectively, of Hj (eq. (Al)) 

panel deflection, positive downward, Z(x,y)e 1WT 
f w \ 2 

v Bxs J 

R 2 

s D xs 

eigenvalue quantities for x- and y-variations, respectively, of assumed beam 
modes (eqs. (A5) and (A10)) 


-i 


- 1 


e p> e q 

V 

X 


width of box, 


w 


B 


’xs 


small difference quantities (see eqs. (A8) and (A14)) 

dummy variable of integration for y £ 

,3 


dynamic pressure parameter, 


2qr 

/3D 


m A 

ratio of panel mass density to air mass density, — — 

PL 

dummy variable of integration for x e 


density of undisturbed airstream 
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length-width ratio of box, 


or 


T 


<Pi 


<pj*(m,n) 


ZB- 


'xs 


wBc 


time 

velocity potential at center of box m,n for mode j, ^(m,n) 
nondimensional part of <pj(m,n) (eqs. (26) and (27)) 


7p velocity potential at center of box m,n due to unit downwash over box 11 , v 

u) frequency of flutter motion 

natural frequencies of panel mode j and of a chosen base or reference 
mode, respectively 

/co bY^ 

S2 = (—-I (1 + ig) flutter eigenvalue parameter 

- 

Q Y 

P 

L J, [ ], { f row, rectangular, and column matrices, respectively 
Subscripts and index numbers: 

B base or reference value 

i, j mode numbers 

m,n streamwise and cross- stream numbered location of a receiving box 

p,q number of half-waves in stream and cross-stream directions, for a panel 

vibration mode 

r,s relative numbered locations of sending and receiving boxes, (eq. (12)) 

te value at trailing edge of panel 

x,y x- and y-directions or variations 

H,v streamwise and cross-stream numbered location of a sending box 


6 



r 


ANALYSIS 


Statement of Problem 

The panel to be analyzed and the coordinate system used are shown in figure 1. 

The panel is a single rectangular panel of length l , width w, and the side edges are 
alined with the remote wind direction. It is a flexible part of an otherwise inflexible sur- 
face that extends at least far enough to the sides to prevent the upper and lower surfaces 
of the panel inducing any aerodynamic effects upon each other. Supersonic flow over the 
upper surface only is considered. The analysis is made on the basis of linear relations 
between small deflections and applied loads. 

Since the primary contribution of the present report is the procedure for obtaining 
the generalized aerodynamic forces for use in a flutter analysis, the governing differen- 
tial equation of motion is written, for simplicity of presentation, for a uniform isotropic 
plate with no in- plane loading 


a 4 z\ a 2 z . / v n 

^ + ^) +m A7? P<X,y ’ T> 


where z is the deflection (positive downward) of the panel, D is the panel flexural 
stiffness, is the mass per unit area of the panel, and Ap is the net perturbation 

pressure (positive downward) arising from the deflection and motion of the panel. Equa- 
tion (1) does not have an explicit term for damping in the structure but any damping con- 


sistent with a linear treatment 
can be introduced. 

For the type of panel deflec- 
tions and aerodynamic forces con- 
sidered, it has not been found 
feasible to make a direct solution 
to the equation of motion. There- 
fore, the often-used procedure is 
followed, of considering the panel 
deflection to be sufficiently well 
approximated by a finite series, 
that is, a linear combination of 
prechosen deflection modes. For 
simple harmonic motion 
z(x,y,r) = Z(x,y)e 1WT , where 



is the circular frequency of Figure 1.- Plan view of panel divided into boxes with coordinate system, 

dimensions, and a forward facing Mach cone with apex at box center 

x m. V 
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vibration and where Z(x,y) can be a complex function to account for phase differences of 
motion over the panel. 


This procedure leads to 

~\ 

z(x,y,T) ~ ^ Hj(x,y, t) 

1 I 

Hj(x,y,r) = q.(r)hj(x,y) = q.e 1Ct,T hj(x,y) 


( 2 ) 


and q^ is the complex amplitude of the generalized coordinate of motion. If the deflec- 
tion mode shapes hj are known that satisfy exactly both the geometric boundary condi- 
tions and the differential equation that corresponds to equation (1) without the aerody- 
namic term, then (the time variation being set aside) 


82 ( H iA) e -icor- 

8T 2 


m A w j 2h j 


( 3 ) 


where on the right-hand side u> takes on its eigenvalue Wj, the normal mode frequency 
associated with hj. The consequences of using shapes hj that fail to satisfy exactly 
the differential equations are examined and discussed later for clamped-edge panels. 
Substitution of equation (2) into equation (1) and application of equation (3) produce the 
result 


^ ^j h j m A C0 j 2 “ X ^j h ) m A w2 " A P( x >y> T ) = 0 (4) 

j j 

in which the frequency go and the combination of q. are unknown, and the perturbation 
pressure Ap still has to be determined. 

Consistent with the use of a series representation of the panel deflection, the pres- 
sure distribution Ap is also used as a series, one term being associated with each 
deflection mode, 

Ap^APj^qjA^j < 5 > 

i i 

where Ap. is separated into the time containing q, and the time invariant Am. 

3 3 3 

The Galerkin method is chosen to form the flutter stability equations. Briefly 
described, the terms of equation (4) are multiplied by a mode shape hq(x,y) and the 
products are integrated over the surface area of the panel. When the modal index num- 
ber i is made successively 1, 2, 3, . . . for all the modes hi employed, the result is 
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a system of equations that requires an equilibrium of work or energy for a condition of 
minimum total (kinetic plus elastic potential) energy of the panel. For modes hj that 
satisfy equation (3) and are orthogonal, each of the equations is of the relatively simple 
form 


^ q.(wj 2 - u) 2 ) n^hjhi dx dy - ^ q. ^Ap.W dx dy = 0 (6) 

j S j S 

and the last term q. \ \ Sp.hj dx dy = Q, . is the generalized aerodynamic force term. 

] J J 3 

s 

In order to obtain the pressures Ap^ = q^ Ap^, the choice is made to work with the 
velocity potential <pj and make use of the relation (compare eq. (3) of ref. 2, for 
example) 


Ap. = -p 


dc? . 
dr 


= -P V 


»(z*) 


x) 8 t) l \9x V J 


(7) 


where p and V are the mass density and velocity of the undisturbed airstream. Each 
of the generalized aerodynamic force terms needed in the dynamic equilibrium 

flutter equations is the integral over the panel surface of work that would be done on 
modal displacement i with unit amplitude due to a deflection in mode j with ampli- 
tude qj; that is, 

Q ij = " ? If hi (“5 + 1 ¥ ^) d(z * )d(w ^ (8) 

S ' ' 


where here and in the following development the subscript index for mode i is not to be 
confused with the unit of imaginaries i = ^-T. The coordinate y ranges from 0 at the 
left edge to 1.0 at the right edge. If the term hi^8<pj/8xj is integrated by parts, the 


result is 


Q ij = ■ xiylo h l( x te'yH x te.y) - ?j(^r - i f 


d(wy) 


(9) 


For elastically supported or free trailing edges, the trailing- edge deflection h^x^ e ,y^ 
could be significant, but for the usual condition of a restrained trailing edge, it is zero, 
and is determined from the second term only as 

Q ij = ' 1 v h ‘) d(a)d(wii) <10) 


Note that the quantity within the parentheses is the complex conjugate of the downwash 
ratio. 
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In most of the foregoing equations, the harmonic time variation e ia,T is 
retained, but except at points where the time derivative 9/ 9 r is indicated, the quantities 
and Q.. can be considered alternatively as their complex amplitudes only. 

In the rest of this section the procedures for computing (1) the aerodynamic influ- 
ence coefficients relating downwash and velocity potential, (2) the velocity potentials, and 
(3) the generalized aerodynamic forces are described. The procedure used for obtaining 
solutions to the system of flutter equilibrium equations is also described. 


Aerodynamic Influence Coefficient 

The panel is considered as divided into a gridwork of equal- size rectangular 
"boxes." The number of boxes in the stream and cross-stream directions are B s and 
B X s> respectively. For reference in the computation procedure, the boxes are numbered 
in sequence beginning at the box nearest the origin of coordinates; in the stream direc- 
tion the index m = 0, 1, 2, . . . Bs - 1, and in the cross- stream direction 
n = 0, 1, 2, . . . Bxg - 1. A second set of box index numbers is needed, and these num- 
bers, designated by p and v, take on values p = 0, 1, 2, . . . B s - 1 and 
v ~ 0, 1, 2, . . . B xs - 1. Thus, reference can be made to the aerodynamic influence 
upon any receiving box m,n due to any sending box p, v. 

The gridwork of boxes is assumed to be sufficiently fine so that the downwash over 
any sending box can be taken as uniformly distributed at any instant, and that the resulting 
pressure perturbation at the center of each receiving box is a sufficiently accurate aver- 
age of the pressure distribution over that box. 

For convenience in the computational procedure, a reference length has been 

chosen as e = , the width of a box. The velocity potential at the center of a 

Bxs 

receiving box m,n due to a uniformly distributed but otherwise unspecified downwash 
w(p,y) over the sending box p,y can be expressed for a simple harmonic time varia- 
tion (compare eq. (34) of ref. 2) as 


<p(m,n) = ew(p,i/) A^(r,s) 


(ID 


where the relative locations fore and aft, and sidewise, respectively, of the two boxes are 
indicated by 


r = m - p | 
s = n - v J 


( 12 ) 
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and from equation (34) of reference 2 


A( p( r ’ s) 'iff 


e -iO(x m -4) cos ^ O 

iM 


[<’ 


m 




*fyn - 


1/2 




(% - «) 2 - - *f 


1/2 


d£ dr? 


(13) 


in which the area of integration v is any part of the sending box \i,v that lies for- 
ward of the forward-facing Mach cone with apex at x m ,y n . (See fig. 1.) Note that the 
integrand is singular along the Mach cone because of the denominator. All boxes and 
parts of boxes lying aft of that Mach cone have a zero contribution to A^(r,s). For 
box m,n the coordinates x m and y n and their respective dummy variables of inte- 
gration | and r\ are based on the reference length e ; that is, x m = j and y n = j. 
The frequency Mach number parameter O is 


- M 2 kg 

n 3 -p- 

p 

where the reduced frequency k e = ^ contains the chosen reference length e . 
As in reference 2, a transformation of variables is made as follows: 



With transformation equations (15), equation (13) becomes 



(14) 


(15) 


(16) 


where the upper and lower limits of the surface integrals, indicated by v^, V 2 > u l> 
and U 2 , are established for any given box by its edges except where the forward Mach 
cone with apex at x m ,y n passes through the box. 

The integral of equation (16) is to be evaluated numerically. The integrand con- 
tains a singularity of order -1/2, at u = ±v, that is, along the Mach lines from x m ,y n . 
In order to minimize the potential inaccuracies in integrating by numerical quadrature 
in the vicinity of that singularity, the numerator of the integrand is replaced by the 
identity 
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cos 


M( u 2 _ v 2) 

2M' 


1/2 


s <cos 


^-( U 2 - v 2 ) 
2M X 


1 / 2 : 


1 ) + 1 


(17) 


Substituting equation (17) into equation (16) and carrying out formally part of the v inte- 
gration gives 


A<p(r,s) 



where the integrand of the remaining v integral is zero at u = ±v and numerical inte- 
gration with good accuracy is possible. All possible forms of the limits uj, U 2 , Vj, 
and v 2 lead to only three forms of equation (18) that are explained with the aid of 
table I. The first form is for the condition V 2 = -Vj = u that can occur only for 
v = n (that is, s = 0), and then only for any portion of a sending box /i, v cut by both 
sides of the Mach cone. For this condition 


L G1 


i p 2 

2j ui 


^(^N fgu du 


(19) 


where Jg is the Bessel function of the first kind of order zero. The second form of 
equation (18) occurs for portions of a box cut by one side of the Mach cone (s * 0), so 

i v 2 

that the limit V£ = u and Vj = 2s - 1 5 1. Since V 2 = u, cos -1 — =0 and 



The third form, which occurs for boxes that are completely ahead of the Mach cone and 
also for portions of boxes ahead of the point where the Mach cone passes out through the 
side of the box, is 


‘OS^re-»k^-e OE -^ + f 


v 2 

vi 


cos 

- v2) 1/2 

- 1 


L2M J 



(u2 - V 2) 


1/2 


dv 


du 


( 21 ) 


The complete A<p(r,s) for any one box fi,v consists entirely of Iq*, Iq 2 > or Iq 3 > or 
a combination of Iqj and or of and I^g. 
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Table I.- Types of integrals and limits of integration for computing A^(r,s) for all 
possible relative locations of box and the Mach cone from box m,n. 





V 1 


0 


0 


0 


2s - 1 


2s - 1 


2s - 1 
2s - 1 


2s - 1 
2s - 1 


2s - 1 


Limits of integration 

v 2 Ul 

1 (2r - 1) e/(3 



1 

0 


1 1 

(2r - 1) a/p 


(2r - 1) a/p 


0 


u (2r - 1) a/p 


2s - 1 


2s + 1 
u 


2s + 1 
2s - 1 


2s + 1 


2s + 1 
(2r - 1) a/p 


2s + 1 (2r - 1) a/p 


A p (r,s) = 0 


u 2 

(2r + 1) a/p 


(2r + 1) a/p 
1 


(2r + 1) a/p 

1 


(2r + 1) a/p 


(2r + 1) a/p 


(2r + 1) a/p 


(2r + 1) a/p 


(2r + 1) a/p 

2s + 1 


(2r + 1) a/p 
2s + 1 


(2r + 1) a/p 


13 





Table I covers all possible relative locations of a sending box and the Mach cone 
from a receiving box. The applicable integral, Iq^> Iq 2’ or *G3’ * or a ^ ox or ^ ox 
segment and the limits of integration v^, Vg, u^, and Ug are listed. The limits 
tabulated are consistent with the reference length e chosen for equation (11) and with 

the transformation of equation (15). The double quantity 2 Iq 3 appears in table I where 

IB 

advantage was taken of symmetry about v 1 = 0. The quantity a = — — — is the length- 

wB s 

width ratio of a box. 

For the unyawed panel with the assumptions that have been made, the aerodynamic 
influence coefficients have right- left symmetry about s = 0; that is, 

A<p(r,s) = A^(r,-s) (22) 

Thus, all values of A^(r,s) that are needed for any box are obtained by considering 
either one of the rear corner boxes as the receiving box. Both the u and v integrals 
are carried out in the computing program according to Gaussian numerical quadrature 
(described in numerous texts, for example, ref. 3) by an available subroutine. For 
sending boxes near the Mach cone from x m ,y n , a five-point integrating rule was found 
to be accurate. For sending boxes remote from the Mach cone, a three-point rule was 
found to be accurate and provided a consequent economy of machine time. 


Velocity Potentials 

Once all possible values of A ( p(r,s) are computed, the total <pj(m,n) at the cen- 
ter of any box m,n for any downwash mode j is a weighted sum of the ^>(m,n) of 
equation (11) particularized to mode j; that is, 


<Pj(m,n) = Ve YY 


UiL-, V 
M v 


(23) 


or, in matrix form 


<p-(m,n) = Ve 


Wj 


A<p(r,s)} 


(24) 


Wi (/i,y) 

The downwash ratios J — are the total time derivatives of Hi 

V J 


w i 1 dHj 8 H) ! 3 H) 5; K , \ 

V V dr 0(Zx) V dr i \dk V 1/ 


(25) 


Of course, a complete array of <pj(m,n) is required to cover all m,n. With a dimen- 
sional constant separated, the complete array is 
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where the nondimensional elements 


n *(m,n) 


<pj*(m,n) are 



( 27 ) 


in which hj is hj^x^,y^. 

The practicality of this box method is largely dependent on the machine time con- 
sumed in the matrix multiplication of equation (27) for multiple modes j. An economical 
way to obtain the matrix of <p j* is described in appendix C. 


Generalized Aerodynamic Forces 

Equation (10) represents a generalized aerodynamic force term as a surface inte- 
gral. With use of the box method, the terms are evaluated by a matrix multiplication. 
With a dimensional constant separated 



w 


BgBxs 


— Q*- 

2 i 


(28) 


the nondimensional part is 



(29) 


To conserve machine core storage, the real and imaginary part of the downwash from 
equation (27) can be used in its complex conjugate in equation (29). 

The complete matrix of elements Qjj is obtained by use of successive rows of the 
downwash conjugate with one row for each mode i and successive columns of < p^* with 
one column for each mode j. 


Solution of the Flutter Determinant 

For the numerical solutions presented herein, it has been assumed that the individ- 
ual modes used in the analysis are orthogonal (that is, with no inertial coupling) and also 
have no significant stiffness coupling between modes. Recent findings, that for clamped- 
edge panels stiffness coupling can be significant, are discussed in the section "Results 
and Discussion." Further, structural hysteretic damping that is characterized by a coef- 
ficient for each vibration mode is assumed. With these assumptions, the system of 



equations that express the dynamic equilibrium of motion based on equation (6) are 


where oq 
mode i 


- an 


i 2 (l + igj) ^ Qij = 0 (i = 1, 2, . . .) (30) 


3 


is the natural frequency of mode i and Mj is the generalized mass of 


p 1 p 1 2 

Mi = Zw J 0 J 0 “AfryJtifcy)] ^ d y 


( 31 ) 


and where represents the distributed panel mass per unit area. In equation (30) 

the quantity gj(«l) is intended to be interpreted as the modal damping coefficient for 
structural hysteretic damping and is made up of two parts 

g x = gi + g 


where g^ can be the assigned or measured values for each mode i of the structure 
and g is a modal-independent mathematical convenience to aid in determining eigen- 
solutions. Let use be made of the asymptotic expression that is applicable for g « 1 
and g i « 1 and exact for either g^ = 0 or g = 0 

1 + ig A + ig ~ (l + ig.)(l + ig) (32) 


(cog/aqw) 2 

If equation (30) is multiplied by — -!■ r, a form convenient for eigensolution results 

Mj(l + igj) 


-£2 + 



1 

1 + igi 


^i 


+ 



a 

(! + igi)kj 2 Mi 


X = 0 


(i = 1, 2, . . .) (33) 


where a>g is any chosen base or reference frequency, 



and 



is the reduced frequency based on the reference length 


l. 


In the set of equations indicated by equation (33), the analyst can readily apply 
values of g i that are not necessarily the same for the various modes i. Of course, 
for the situation where the g^ for all modes are equal, g = gg for illustration, flutter 
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characteristics that are identical in every respect are computed by either of two pro- 
cesses: (1) Assign g^ = gg for all modes and determine the crossing point g = 0; or 
(2) Assign gj = 0 for all modes and determine the crossing point g = gQ. 

For panels with a uniform mass distribution m^, the quantity a/Mi in equa- 
tion (33) can be replaced as follows: 


— _ 1 
M i ^ 


(34) 


where 


2 

l B s Bx S 

Ml * ‘ $0 So h ‘ 2,s } 

1 _ Pi 
M m A 


(35) 


the ratio of air mass to panel mass being 1/p. 

As is commonly the case with the unsteady air forces, the flutter boundary cannot 
be computed directly, and an indirect method of computing and cross-plotting is neces- 
sary. In advance of a flutter solution for any given panel and Mach number, choices are 
made for B S ,B XS and for the number of modes in the analysis. The downwash quantities 
9h^/9x and lq at the center of each box and the modal frequencies oq are computed 
or obtained from experiment for each mode i and arranged systematically for use in 
the matrix multiplications. Values of g^ are assigned and the quantities Mj and a 
^or Mj*, 1/ j u., and a * for panels with uniform m A j are computed, except that pro- 
vision for varying the air density in a. (or 1/ p) is retained. 

After a choice of the reduced frequency based on experience is made, the 
matrix of values A^(r,s), the (£>j(m,n) at each box center for each mode, and the gen- 
eralized aerodynamic forces are computed and employed in the flutter determinant, 

that is, the determinant of the matrix of the coefficients of qj from equation (33). For 
each of a number of values of the air density, the set of complex eigenvalues fi, the 
eigenfrequencies oj, the values of g, and the associated stiffness parameters 

a) 3 L j a) ^ 

= — are computed. There is a value of each quantity for each chosen panel 

mode. The normalized eigenvectors for each eigenfrequency can also be obtained, and 
can be used to determine the proper plotting of curves through the computed results for 
a range of or of a*. By plotting curves of g against 1/ p. for a sufficient range 
of air densities, the existence or nonexistence of an eigenvalue with g = 0 for some air 
density is established for each mode. An associated curve of the stiffness parameter 
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a>l ijV or j a s is plotted against mass ratio - one curve for each eigenvalue. 

(The use of the speed of sound a s rather than V permits direct comparison of flut- 
ter boundaries for various Mach numbers.) 

On a cross plot of 1/ 1 ± against stiffness parameter, points can now be placed to 
represent g = 0. Each point is on a separate stability boundary. 

By successive choices of based on the results of previous choices, a series of 
points can be placed on the plot of 1/ ji, against stiffness parameter, and stability bound- 
aries can be drawn through the points. Four different types of stability boundaries have 
been found and three types are illustrated in figure 2. For each type, a boundary for 
= 0 and one for g^ = 0,01 (as illustrative of some small value) are shown, and the 
arrows point in the direction of increasing reduced frequency kj. The boundary of fig- 
ure 2(a) predicts that if g^ were zero, a fixed panel thickness would be needed as the 
air density tends toward zero, but if g^ were not zero, panel thickness would go toward 
zero as air density does. 


For the boundary of figure 2(b), the unstable 

region shrinks with increasing g i and vanishes for 

some value that may be extremely small or as large 

as g i = 0.05 or possibly larger, depending upon the 

panel and flow parameters. The crossing point of 

the boundaries for g. = 0 in both figures 2(a) and 
1 1 

2(b) through - = 0 occurs at a value of k^ for 
which the imaginary part of the Qj^kj) passes 
through zero and the panel motion has a pure single 
degree of freedom. For the boundary of figure 2(c), 
the flutter motion is strongly coupled. Small 
increases of gj have little or no effect, and this 
small effect can be either stabilizing or destabil- 
izing, depending on the particular panel and flow 
parameters. An important point regarding the 
boundary of figure 2(c), as well as the dashed bound- 
ary of figure 2(a), concerns their resemblance to a 
parabola. For any boundary or portion of a bound- 
ary which can be represented by a parabola 

described by -i- = 

density and the panel thickness ratio t /l to pre- 
vent flutter are related by the formula 


o> i 


— i- times a constant, the air 



= Constant. 


Such a relationship is contained 




Figure 2.- Three general types of stability 
boundaries for g, = 0 and for gj = 0.01. 
Arrows indicate direction of increasing 
reduced frequency kj. 



in the panel flutter parameter (Ej8/q)^(t/Z), which has been evolved and used by sev- 
eral investigators. The fourth possible type of boundary found from a flutter determinant 
is one that falls entirely in the negative 1/ p. region, the positive 1/ p region being 
stable with respect to it. 

For any given panel and flow parameters, a set of stability boundaries is computed 
according to the number of modes used in the analysis. The flutter boundary separates 
the region that is stable with respect to all stability boundaries from the region that is 
unstable with respect to at least one stability boundary. Care should be taken to estab- 
lish that the flutter boundary is converged; that is, enough modes have been used in the 
analysis so that additional modes do not alter the flutter boundary in any important 
respect. By the Galerkin method alone, it is not usually possible to establish substantial 
convergence with mathematical rigor. When the incremental effect of additional modes 
is not only small relative to the magnitude of the calculated quantities, but is also 
decreasing at a substantial rate as modes are added, near convergence is commonly 
assumed. 


RESULTS AND DISCUSSION 

Results that have been calculated by means of the present analysis are presented in 
a number of figures to scales large enough to permit the reading of values from the flut- 
ter boundaries for design purposes. From past work it was established, at least for 
wide panels (Z/w less than about 1.0), that for the low supersonic speed range (M less 
than about i/2), approximate aerodynamic forces, such as those of piston theory and static 
strip theory (the latter is also known as the Ackeret value), do not give valid analytical 
results. Nevertheless, because of their ease of applications, and because of their applica- 
tion to the higher supersonic speed range, interest in the use of approximate aerodynamic 
forces has continued. This continuing interest has led to the as-yet- limited finding that 
for panels that are not wide (Z/ w greater than about 1.0), approximate aerodynamic 
forces result in flutter boundaries that agree at least fairly well with the boundaries from 
the present analysis even in the low supersonic speed range. (See fig. 2 of ref. 4.) On 
the figures that apply to ^ ^ 1.0, therefore, the flutter boundary as obtained from the 
roots of the closed-form expressions of reference 5 for simply supported panels and of 
reference 6 for clamped- edge panels (with minor changes for the latter as explained sub- 
sequently) is also shown. Reference 5 demonstrates that with static strip-theory aero- 
dynamics on the simply supported panel, the critical flutter mode from the closed-form 
solution involves a single half- sine wave as its cross -stream variation. The closed-form 
expressions of reference 6 for the clamped- edge panel were developed on the basis of an 
assumed single half-wave in the cross- stream direction. Since the natural frequencies 
for multiple half-wave modes are always higher than those for the single half-wave modes, 
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and since the aerodynamic cross coupling between single and multiple half-wave modes 
has been observed to be at least fairly small in comparison with the direct coupling 
terms, all the results presented in this report are for single half-wave cross- stream 
modes. 

Figures 3 and 5 present the flutter boundaries in plots of the mass ratio l/(i as 
the vertical coordinate and the stiffness parameter V as the horizontal coordinate. 

Flutter boundaries are presented for M = 1.3 for simply supported panels with length- 
width ratios of 0, 1/4, 1/2, 1, 2, 4, and 10, and for clamped-edge panels with the same 
length-width ratios except 10. All the results were obtained or spot-checked with 
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(a) Pinned edges; ^ = 0 (two-dimensional case). 

Figure 3.- Flutter boundaries and other noncritical stability boundaries for plane, unyawed, unstressed, isotropic, 

rectangular panels with M = 1.3. 


20 



B s = 40 and B xs = 10, except that Bx S = 1 for ~ = 0. Values of the structural 
damping coefficient are equal for all modes. 

Appendix A gives the expressions for the natural modes from which the values of 
downwash ratio for each box center and the modal frequencies were computed. 

Panels With Simply Supported Edges 

Figures 3(a) to 3(g) present the results for the simply supported panels. Fig- 
ure 3(a) contains portions of the flutter boundary and of three noncritical stability bound- 
aries for the two-dimensional panel ^ = oj . The boundaries shown are parts of more 
complete curves of the type shown in figure 3 of reference 7. The solid curves are the 
boundaries for g i = 0, and the short-dashed curves are for g i = 0.01. The unstable side 
of each boundary is to the left. Consequently, the boundary farthest to the right forms 
the critical or flutter boundary, and is indicated by tick marks on the solid and short- 
dashed curves. The predominant natural-mode components for the four stability 
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Figure 3.- Continued. 
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boundaries are given by the numerals 1 to 4. The mode-1 boundary is the critical one for 
this panel. For the three points labeled by the letters A, B, and C, additional information 
is given in the table of the figure, including the parameters 1/ p, g^ kj, co^Z/V, 
w/wj, and the six eigenvector components of the flutter mode. The eigenvectors with 
their small imaginary parts show that the flutter motion is nearly a standing wave with 
the shape of the first vacuum mode; that is, the flutter motion is nearly the single degree 
of freedom described by the first mode. 

On figure 3(a) for two-dimensional panels, in addition to the flutter and noncritical 
stability boundaries, there are four hyperbolas (long-dashed lines) obtained from the U.S. 
Standard Atmosphere (ref. 8) that apply to aluminum panels at four altitudes, namely, sea 
level, 10 000, 25 000, and 50 000 feet. The hyperbolas for other altitudes and other panel 
materials can be determined as desired. The intersection of the appropriate hyperbola 
with the flutter boundary determines the thickness ratio t/l required to prevent flutter 
for a given panel material and altitude. For simply supported aluminum panels at sea 
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(c) Pinned edges; ^ 
Figure 3.- Continued. 
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level, the boundary for gj = 0 specifies a minimum thickness ratio of = 0.0181, 

f 

whereas the gj = 0.01 boundary specifies — = 0.0146. For aluminum panels at 

If 

50 000 feet on the gi = 0.01 boundary, y- = 0.0087. 

Figure 3(b) gives the results for simply supported panels with = -j. The flutter 
boundary is predominantly mode 1, that is, nearly a standing wave (see table on fig. 3(b) 
for points A, B, and C), and the noncritical stability boundaries included are predom- 
inantly modes 2, 3, and 4. For aluminum panels at sea level with g. = 0, ^- = 0.0138; 
this value is a 24-percent reduction from the corresponding value for the two- 
dimensional panel. 

Figure 3(c) gives the results for panels with = ~. For the 1/ /i range of the 
figure, the mode- 2 boundary is now the flutter boundary for gj = 0 and also for 



(d) Pinned edges; y = 1. 
Figures.- Continued. 
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g. = 0.01 except in the very lowest range of l//x. For aluminum panels with g. = 0 at 

f 

sea level, — = 0.0095 and is virtually a constant value for all higher altitudes. With 
g. = 0.01, however, y = 0.0082 at sea level and decreases to y = 0.0059 at 50 000 feet. 

it i 

The three points labeled A, B, and C correlate with the quantities given in the table of 
figure 3(c), and the tabulated values show that the flutter motion is nearly a single- 
degree- of-freedom standing wave in the second natural mode. Point B corresponding to 
the higher air density contains the larger components of modes 1 and 3, and these com- 
ponents are so phased that the aft portion of the panel has a slightly larger amplitude 
than the forward half of the panel. The mode-1 stability boundary, labeled with the 
numeral 1, is no longer the critical boundary unless g i should become as large as about 
0.03. In figure 3(c) the character of the mode-1 boundary for g i = 0 (the boundary for j 
g^ = 0.01 is virtually coincident) is entirely different from those of figures 3(a) and 3(b). 
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(e) Pinned edges; ^ ■ 2. 
Figure 3.- Continued. 
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Instead of intersecting the horizontal axis, it departs from the origin in the manner of a 
near parabola. 

Figure 3(d) gives the results for panels with ^ = 1. The mode- 3 boundary now 
forms the flutter boundary for g^^ = 0 and g i = 0.01, and the flutter mode is predomi- 
nantly the third natural mode (see points A, B, and C and the table), small proportions of 
modes 2 and 4 being present. If the structural damping coefficient could somehow be as 
large as about 0.03, the mode-1 boundary would become the flutter boundary since it is of 
the type that is only slightly affected by an increase of gj to 0.03. The flutter motion 
for the mode-1 boundary is characterized by having the natural mode 1 predominate, 
having a large proportion of mode 2 present, and having the proportions of the higher 
modes decrease. Refer to point D on figure 3(d) and its column in the table. The motion 
is nearly a standing wave, as evidenced by the small imaginary parts of the generalized 
coordinates, and the aft part of the panel has a substantially larger amplitude than the 
forward part. For aluminum panels at sea level with g^ = 0, jf- = 0.0060, and for 

g i = 0.01, ~ = 0.0056. The parabolic boundary obtained from equation (9) of reference 5 
is included as a matter of interest and it is observed to fall close to the flutter boundary 
of the present analysis for g i ~ 0.03. 

Figure 3(e) gives the results for panels with ^ = 2. With g i = 0 the mode- 5 
boundary is critical for most of the range of 1/ p (less than about 0.13), and the mode-1 
boundary is critical for higher values of 1/ p. With g i = 0.01, the numbers labeling the 
segments of short-dashed curves show that the mode-1 boundary is critical for the high 
and low ranges of 1/p. In the intermediate range of l/p, the flutter boundary is 
formed partly by the mode-6 boundary and partly by the mode- 5 boundary. The asso- 
ciated tabulated values show that at point A on the mode-1 boundary, the flutter motion 
is strongly coupled; natural modes 1 and 2 are present in nearly equal parts, natural 
mode 3 is present in a significant amount, and the other modes are present in small 
amounts. The mode-1 boundary of figure 3(e) has this characteristic throughout the 
range of 1/ p of the figure. At points B and C on the mode- 5 boundary and point D on 
the mode-6 boundary, the motion is dominated by the respective natural mode, moderate 
to small amounts of other modes being present. For aluminum panels at sea level with 
g^ = 0, j = 0.0035 (mode- 5 boundary) and with g^ = 0.01, ~ = 0.0034 (mode-1 bound- 
ary). The parabolic boundary from equation (9) of reference 5 is so closely coincident to 
the mode-1 boundary that it is not drawn. 

Figure 3(f) displays the results for panels with ^ = 4. With g i = 0, the flutter 
boundary is formed, beginning at the lowest value of 1/p, by segments of the mode- 8, the 
mode-9, and the mode-2 boundaries. With g^ = 0.01, the flutter boundary is made up of 
segments of the mode-2, the mode-12, the mode-11, the mode- 10, and the mode-2 bound- 
aries. On the mode- 2 boundary, as indicated by the tabulated values for point A, the 
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flutter motion contains large components of several of the natural modes, the mode- 2 
component being the largest. The imaginary parts of the components indicate that the 
flutter motion is not (even approximately) a standing wave. On the other boundary seg- 
ments labeled 8, 9, 10, 11, and 12, those respective natural modes are predominant in the 
flutter motion with moderate to small proportions of other modes being present for 
points B, C, D, and E as indicated in the table of figure. 3(f). For aluminum panels at sea 
level, j = 0.0020 for both gj = 0 and gj = 0.01, since the mode- 2 boundary is not sig- 
nificantly affected by structural damping. In figure 3(f) the parabolic boundary from 
equation (9) of reference 5 falls near the mode-2 boundary and on the conservative side 
and calls for a greater thickness to prevent flutter 

Figure 3(g) shows the results for panels with ^ - 10. The flutter boundary is 
formed by the mode- 7 boundary throughout the range of 1/ p shown. For aluminum 
panels at sea level, the flutter boundary corresponds to ^- = 0.00097 for g; = 0 and to 
- = 0.00096 for g^ = 0.01. The mode-9 boundary (natural modes 3 to 14 used in the anal- 
ysis) and the mode- 11 boundary (natural modes 5 to 16 used in the analysis) are in the 
unstable region very close to the flutter boundary. Stability boundaries associated with 
modes 2, 4, and 6 (not shown) also lie within the unstable region, the mode- 2 boundary 
being the farthest from the flutter boundary. For the three points A, B, and C on the flut- 
ter boundary, parameters and generalized coordinates are given in the table of figure 3(g). 
The value of is markedly larger than that for any other boundary in figure 3 (except, 

of course, the mode- 12 and mode- 11 boundaries of figure 3(f). This large value of q^ 
leads naturally to the question of whether 12 modes are sufficient in this particular modal 
solution to provide a substantially converged result. Therefore, spot checks were made 
with 18 modes. Since more than 12 modes could be used only outside the main program 
and the required intermediate work was extensive, spot checks were made for only two 
values of reduced frequency. With 18 modes the flutter boundary moved inward to the 
left by an amount that would reduce the thickness required to prevent flutter by about 
8 percent. The mode- 7 boundary still formed the flutter boundary, and other stability 
boundaries also moved and maintained their relative locations. At the two widely sepa- 
rated points obtained on the mode- 7 flutter boundary (see the last 2 columns in table of 
fig. 3(g)), the generalized coordinates decreased monotonically for the modes above q^. 
For the 18th mode, jq-^gj ~ 0.02, whereas for mode 12, (qx2| ~ 0-14. The result is, there- 
fore, thought to be substantially converged, although rigorous proof is lacking. The 
parabola of long dashes from equation (9) of reference 5 corresponding to a value of 
X = 51 390 is close to the flutter boundary of the present analysis. The associated flut- 
ter frequency ratio — -1.18 and falls just above the fourth natural mode -j- - 1.15. 

Additional information for this panel is given by figure 4 that presents the flutter 
mode shapes as obtained from the present analysis and from an automated solution of 
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(f) Pinned edges; ^ = 4. 
Figure 3.- Continued. 
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(0.097 - O.lOi) x 10-1 

-0.10 + 0.036i 

1.0 

(-0.50 + 0.28i) x 10-2 

(-0.13 + 0.042i) x lO’l 

(0.61 - 0.73i) x io-2 

-0.096 + 0.023i 


(0.21 - 0.24i) x 10-2 

(-0.99 + 0.331) x IQ" 2 

(0.55 - 0.44i) x io-2 
(-0.89 + 0.20i) x io-2 
(0.13 - O.lli) x iq-2 


Point A 

Point B 

Point C 

* 

* 

0.464 

0.302 

0.174 

0.662 

0.057 

0 

0 

0 

0 

0 

4.5 

3.8 

3.0 

4.5 

1.5 

3.415 

2.825 

2.20 

3.52 

1.13 

1.3177 

1.344 

1.37 

1.28 

1.33 

0.052 - 0.000 15i 

0.047 - O.OOOSSi 

0.043 - 0.00079i 

0.10 + 0.015i 

0.073 + 0.0020) 

-0.111 - 0.0022i 

-0.099 + 0.00026i 

-0.092 + 0.00078i 

-0.21 - 0.033i 

-0.16 - 0.0046i 

0.19 + 0.0078i 

0.17 + 0.0022i 

0.16 + 0.00066i 

0.33 + 0.057i 

0.26 + 0.0088i 

-0.30 - 0.027i 

-0.28 - 0.014i 

-0.26 - 0.0078i 

-0.48 - 0.089i 

-0.40 - 0.016/ 

0.48 + 0.062i 

0.45 + 0.038i 

0.42 + 0.024i 

0.67 + 0.12i 

0.60 + 0.026i 

-0.74 - O.lli 

-0.73 - 0.080i 

-0.71 + 0.0571 

-0.87 - 0.12i 

-0.85 - 0.029i 

1.0 

1.0 

1.0 

1.0 

1.0 

-0.81 + 0.23i 

-0.83 + O.I9i 

-0.85 + 0.15i 

-0.88 + 0.17i 

-0.85 + 0.041/ 

0.44 - 0.21i 

0.45 - 0.17i 

0.47 - 0.13i 

0.60 - 0.21i 

0.55 - 0.043i 

-0.25 + 0.13i 

-0.26 + O.lOi 

-0.26 + 0.076i 

-0.38 + 0.16i 

-0.33 + 0.028/ 

0.15 - 0.074i 

0.14 - 0.057i 

0.15 - 0.041i 

0.24 - O.lli 

0.20 - 0.017i 

-0.11 + 0.058i 

-0.11 + 0.043i 

-0.11 + 0.0311 

-0.15 + 0.070i 

-0.13 + O.OlOi 




0.10 - 0.0451 

0.083 - 0.0063i 




-0.071 + 0.031i 

-0.057 + 0.0043i 




0.050 - 0.021i 

0.040 - 0.00251 




-0.037 + 0.015i 

-0.029 + 0.0037i 




0.027 - O.OlOi 

0.021 + 0.00036i 




-0.022 + 0.0087/ 

-0.017 - 0.00026) 


Discussed in text only. 
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(g) Pinned edges; ^ - 10. 
Figure 3.- Concluded. 
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Figure 4.- Flutter mode shapes from present analysis with 12 and with 18 modes, M = 1.3, and 
= 4.5, and from equation (9) of reference 5 for a simply supported rectangular panel with 
1 / w = 10, g,- = 0. 

equation (9) of reference 5. The latter result is real only and is associated with the point 
at a length-width ratio of 10 on the solid curve of figure 2 of reference 4. From the pres- 
ent analysis, both the real (inphase) and imaginary (quadrature) parts as obtained from 
using 12 modes and also from using 18 modes are shown. The two results are for 
kjr = 4.5. The 12-mode result corresponds to point A in figure 3(g) and the first column 
of the table. The 18-mode result corresponds to the tabulated column second from the 
right. An examination of those two columns and the corresponding pairs of curves of 
figure 4 reveals their moderate differences. The closed-form solution based on refer- 
ence 5 is like the solutions from the present analysis in that the point of maximum deflec- 
tion is at about 93 percent of the panel length, and very little motion occurs over the for- 
ward half of the panel; the solution differs in that the signs forward of the maximum 
deflection point do not alternate. A Fourier analysis, with 101 points and 0 ^ n ^ 50, of 
the closed-form deflection curve has given the results of table II. The largest proportion 
present is for p = 3 half-waves, a gradual and monotonic decrease of proportions for p 
increasing above 3; and there is an unfailing alternation of sign for even and odd p. The 
proportions of sine waves in table II can be compared with the generalized coordinates 
in the columns of the table in figure 3(g). 


Panels With Clamped Edges 

Figures 5(a) to 5(f) present the results for clamped- edge panels. In each of the six 
figures the solid curve with tick marks attached is the flutter boundary for g. = 0, and 


29 


TABLE II.- PROPORTIONS OF SINE WAVES WITH n HALF-WAVES PRESENT IN 
THE DEFLECTION SHAPE LABELED "REF. 5" IN FIGURE 4 


Number of half-waves, 
n 

Proportions present 

0 

0 

1 

.544 

2 

-.901 

3 

1.000 

4 

-.905 

5 

.722 

6 

-.536 

7 

.382 

8 

-.267 

9 

.186 

10 

-.130 

11 

.0923 

12 

-.0660 

13 

.0478 

14 

-.0351 

15 

.0261 

16 

-.0196 

17 

.0149 

18 

-.0115 

19 

.00895 

• 

49 

.0000580 

50 

-.0000507 


the curve of short dashes with tick marks is the flutter boundary for g. = 0.01. The flut- 
ter region is to the left of the flutter boundary in all the figures. In addition to the flutter 
boundary, some of the noncritical stability boundaries that also result from the eigensolu- 
tion are shown without tick marks so that the progressive alteration of the flutter boundary 
with increasing values of 1 / w can be more readily appreciated. Also given in each fig- 
ure are the hyperbolas that apply to aluminum panels at sea level and at 10 000, 25 000, 

and 50 000 feet in a standard atmosphere (ref. 8). In figures 5(d) to 5(f) with ~ g 1.0, the 

w 

closed-form solutions based on the Ackeret value of static aerodynamic forces are also 
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presented as the labeled parabolas. The solutions were obtained essentially as described 
in reference 6 except for the minor differences that (1) the cross-stream variation of 
mode shape was assumed to be the natural vibration mode of a uniform clamped-end 
beam with a single half-wave instead of a 1.0-minus-cosine type of half-wave (see 
eq. (8.37) of ref. 6), and (2) the dynamic pressure parameter X was defined as in ref- 
erence 5 as X = and is consistent with the assumed aerodynamic forces. 


0D 


Figure 5(a) applies to ^ = 0 (the two-dimensional case) and the solid lines shown 
are essentially identical to portions of those from figure 3 of reference 7. The predom- 
inant modes are indicated by the number for each boundary. For the three points labeled 
by the letters A, B, and C on the flutter boundaries, information is given in the table. 
Clearly, the flutter mode is virtually a standing wave with a small amount of the second 
natural mode present along with the predominant first natural mode. For aluminum 
panels at sea level, the thickness ratio for gj = 0 is |- = 0.0098, and for gj = 0.01, 


Quantity 

Point A 

Point B 

Point C 

1/M 

0.0441 

0.0527 

0.01087 

9i 

0 

0.010 

0.010 

k ( 

0.745 

0.66 

0.40 

I Pjl/V 

0.7583 

0.6690 

0.4007 


0.982 

0.987 

0.998 

h 

1.0 

1.0 

1.0 

f)2 

-0.050 + 0.0046i 

-0.077 + 0.0058i 

-0.042 + 0.00121 

d3 

I0.13 - 0.093i) x 10‘ 2 

(0.28 - O.lli) x 10-2 

(0.89 - 0.17i) x 10-3 

04 

(-0.12 + 0.024i) x 10-2 

(-0.18 + 0.030i) x IQ' 2 

(-0.10 + 0.005i) x 10-2 

05 

(0.90 - 0.91i) x 10-4 

(0.18 - O.lli) x 10-3 

(0.61 - 0.15i) x io-2 

% 

(-0.13 + 0.037i) x io-3 

(-0.21 + 0.045i) x 10"3 

(-0.11 + 0.006i) x iq-3 






Constant values of for 

aluminum panels at four altitudes 



V 


(a) Clamped edges; ^ = 0 (two-dimensional case). 

Figure 5.- Flutter boundaries and other noncritical stability boundaries for plane, unyawed, unstressed, isotropic, 

rectangular panels with M = 1.3. 
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j- = 0.0086. For aluminum panels at 50 000 feet on the gj = 0.01 flutter boundary, 
t = 0.0058. 

Figure 5(b) applies to ^ When compared with figure 5(a), the unstable region 
has shrunk, the mode-1 boundary having moved to the left considerably more than the 
mode- 2 boundary has. The result is that in the lower range of 1/ \i, the mode- 2 bound- 
ary forms the flutter boundary for both = 0 and g i = 0.01. As indicated in the table 
of figure 5(b), points A and B are predominantly mode 1, some mode 2 being present; and 
point C is predominantly mode 2, small amounts of modes 1 and 3 being present. For 
aluminum panels at sea level, j = 0.0081 for g A = 0, and j = 0.0070 for gj = 0.01. 

Figure 5(c) applies to ^ For g^ = 0 the mode- 2 boundary forms the flutter 
boundary throughout the 1/ / u range shown. For g. = 0.01, the mode-2 boundary is 
critical for about 0.01 < 1 / \x < 0.10. The parameters listed for points A, B, C, and D 
show that the flutter mode is predominantly natural mode 2 and has small to moderate 
amounts of modes 1 and 3. The amounts of modes 1 and 3 that are present increase with 
increasing values of l//r. The mode-1 boundary has changed its character from that for 


Quantity 

Point A 

Point B 

Point C 

1/p 

9i 

H 

Wil/M 

u/ui 

9l 

q 2 

q 3 

94 

95 

96 

0.0566 

0 

0.635 

0.6475 

0.981 

1.0 

-0.091 + 0.0070i 
(0.37 - 0.16i) x 10“2 
(-0.22 + 0.040i) x HP 
(0.23 - 0.14i) x io-3 
(-0.24 + 0.061i) x 10-3 

0.06173 

0.010 

0.575 

0.5790 

0.993 

1.0 

-0.12 + 0.0086i 
(0.64 - 0.20i) x 10-2 
(-0.31 + 0.050i) x 10-2 
(0.41 - 0.18i) x 10-3 
(-0.35 + 0.0731) x 10’3 

0.0192 

0.010 

1.2 

0.4447 

2.70 

-0.063 + 0.018i 
1.0 

-0.036 + 0.0021i 
(0.12 - 0.058i) x 10' 2 
(-0.15 + 0.023i) x 10-2 
(0.15 - 0.087 i) x 10-3 



v 

(b) Clamped panel; ^ 
Figure 5.- Continued. 
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— = ~ and 0; in figure 5(c) it is approximately a parabola with its apex at the origin. 

W 4 

Furthermore, the flutter mode is a coupled one, mainly of modes 1 and 2, formed by a 
near coalescence of those two flow-affected frequencies. If gi should be greater than 
0.03, the mode-1 boundary is the flutter boundary for the range of 1/ j u. of the figure. 

For aluminum panels at sea level, j- = 0.0068 for gi = 0, and j = 0.0059 for gi = 0.01. 

Figure 5(d) applies to — = 1. For g^ = 0 the mode-3 boundary forms the flutter 

w 

boundary in the 1/p. range shown. For gi = 0.01, the flutter boundary is formed by the 
mode-1, -3, and -4 stability boundaries. For gi greater than about 0,03, the mode-1 
boundary is the flutter boundary. For the five points labeled A to E on the flutter bound- 
aries, the proportions of the first six natural modes are given in the table of figure 5(d). 


Ouantity 

Point A 

Point B 

Point C 

Point D 

Point E 

1/M 

0.0725 

0.0775 

0.0370 

0.0120 

0.0828 

9i 

0 

0.010 

0.010 

0.010 

0 

H 

1.46 

1.25 

1.2 

0.8 

0.5 

wjt/V 

0.5833 

0.5022 

0.471 

0.3101 

0.4651 

C0/C0] 

2.50 

2.48 

2.55 

2.58 

1.075 

9l 

-0.14 + 0.050i 

-0.22 + 0.059i 

-0.12 + 0.0341 

-0.094 + 0.016i 

1.0 

q 2 

1.0 

1.0 

1.0 

1.0 

-0.30 + 0.018i 

93 

-0.089 + 0.015i 

-0.12 + 0.014i 

-0.069 + 0.0073i 

-0.049 + 0.0017i 

0.034 - 0.0060i 

94 

(0.49 - 0.440 x io‘2 

(0.95 - 0.47i) x 10-2 

(0.35 - 0.20i) x 10-2 

(0.19 - 0.0420 x 10‘2 

(-0.90 + 0.15i) x 10‘2 

q 5 

(-0.36 + 0.14i) x 10-2 

(-0.55 + 0.14i) x io-2 

(-0.30 + 0.62i) x 10-2 

(-0.21 + 0.014i) x 10-2 

(0.21 - 0.058i) x 10-2 

96 

(0.47 - 0.58 i) x iq-3 

(0.96 - 0.64i) x iq-3 

(0.38 - 0.28i> x 10-3 

(0.21 - 0.058i) x iq-3 

(-0.10 + 0.023i) x iq-2 



(c) Clamped panel; ^ = 


Figure 5.- Continued. 
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The mode-1 boundary is characterized by a strongly coupled motion. (See point A.) 
The other boundaries each have a predominant mode, moderate to small proportions of 
adjacent modes being present. For aluminum panels at sea level, = 0.0048 for 

t V 

gi = 0, and — = 0.0043 for gj = 0.01. The parabolic boundary from the closed-form 

0 

solution is included and falls on the unconvervative side of the mode-1 boundary. 

Figure 5(e) applies to J- = 2. For g. = 0, the flutter boundary is formed by the 
mode-1 and the mode- 5 boundaries, the mode-6 and mode-7 boundaries being near the 


Quantity 

Point A 

Point B 

Point C 

Point D 

Point E 

1/M 

0.0966 

0.1104 

0.0940 

0.0436 

0.0163 

9| 

0 

0.010 

0.010 

0.010 

0.010 

kf 

2.12 

0.628 

2.96 

1.8 

2.1 

wjf/V 

0.6050 

0.5188 

0.5249 

0.502 

0.3606 

U/U 1 

3.50 

1.21 

5.63 

3.59 

5.82 

h 

0.032 - 0.018i 

1.0 

-0.013 + 0.0025i 

0.019 - 0.0063i 

(-0.34 - 0.073i) x io' 2 

92 

-0.19 + 0.082i 

-0.64 + 0.068i 

0.029 - 0.016i 

-0.13 + 0.049i 

(0.73 + 0.36i) x io-2 

93 

1.0 

0.16 - 0.034i 

-0.15 + 0.087 i 

1.0 

-0.068 + 0.016i 

94 

-0.14 + 0.032i 

-0.032 + 0.0090i 

1.0 

-0.091 + 0.014i 

1.0 

95 

0.011 - O.Olli 

(0.99 - 0.32i) x 10-2 

-0.13 + 0.035i 

(0.58 - 0.42i) x iq-2 

-Q.044 + 0.0026i 

96 

(-0.74 + 0.32i) x 10' 2 

(-0.34 + 0.12i) x 10'2 

0.0010 - 0.012i 

- 

(-0.50 + 0.13i) x 10-2 

(0.21 - 0.071D x io-2 


Constant values for 

aluminum panels at four altitudes 



(d) Clamped panel; ~ - 1. 
Figure 5.- Continued. 
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mode- 5 boundary. For g^ = 0.01, the flutter boundary is formed by portions of the 
mode-1, -6, -7, and -8 stability boundaries. (The mode-7 and -8 boundaries are virtually 
coincident over part of the range.) For g i as large as about 0.03, the mode-1 boundary 
is the flutter boundary. The parameters listed included the eigenvectors for point A 
(coupled motion) and also for points B, C, D, and E where one mode is predominant for 
each. For aluminum panels at sea level, = 0.0031 for both g^ = 0 and 0.01. The 
parabolic boundary from the closed-form solution falls close to the mode-1 boundary. 

Figure 5(f) applies to ^ = 4. For g^ = 0, the flutter boundary is formed by por- 
tions of the mode-2, the mode- 10, and the mode-9 stability boundaries. The mode-11 and 
mode- 12 boundaries (and other higher mode boundaries not shown) lie to the left of the 


Quantity 

Point A 

Point B 

Point C 

Point D 

Point E 

1 /m 

0.1398 

0.0753 

0.0400 

0.0280 

0.0274 

9i 

0 

0 

0.010 

0.010 

0.010 

k l 

1.12 

3.4 

3.6 

4.8 

3.6 

wg/V 

1.063 

0.9829 

0.768 

0.624 

0.590 

UJ/Cdi 

1.053 

3.46 

4.69 

7.69 

6.10 

9l 

1.0 

(0.38 - 0.33i) x 10-2 

(-0.17 - 0.0460 x lo- 2 

(-0.43 - 0.130 x 10 3 

(0.45 - 0.420 x 10-3 

q 2 

-0,88 + 0.24i 

-0.015 + 0.0036i 

(0.30 - 0.13i) x 10-2 

(0.53 - 0.250 x 10-3 

(-2.1 - 0.150 x io-3 

q 3 

0.29 - 0.15i 

0.027 - 0.014i 

(-1.0 + 0.200 x 10-2 

(-1.9 + O.lOi) x 10-3 

(3.2 - 0.380 x 10-3 

q 4 

-0.063 + 0.040i 

-0.16+ 0.081i 

(1.5 + 0.0540 x 10-2 

(2.3 - 0.0590 x io-3 

(-0.91 + 0.370 x io-2 

% 

0.019 - 0.013i 

1.0 

-0.12 + 0.044i 

(-0.65 + 0.220 x 10-2 

(0.71 + 0.410 x 10-2 

96 

(-0.62 + 0.49i) x 10' 2 

-0.13 + 0.042i 

1.0 

(0.74 + 0.460 x 10-2 

-0.11 + 0.022i 

97 


(0.91 - 1.2i) x 10-2 

-0.091 + 0.019i 

-0.082 + 0.024i 

1.0 

98 


(-0.94 + 0.390 x io-2 

(0.54 - 0.470 x 10"2 

1.0 

-0.080 + 0.0098i 

99 


(0.15 - 0.230 x 10“2 

(-0.69 + 0,170 x 10-2 

-0.063 + 0. 0095 i 

(0.48 - 0.260 x io-2 

9 10 


(-0.21 + O.lOi) x 10-2 

(1.0 - 0.990 x 10-3 

(0.33 - 0.230 x 10-2 

(-0.65 + 0.960 x 10-2 

911 


(0.43 - 0.740 x 10-3 

(-1.7 + 0.480 x 10’ 3 

(-0.52 + 0.0860 x io-2 

(0.97 - 0.57 0 x 10-3 

912 


(-0.68 + 0.410 x iq-3 

(0.35 - 0.340 x ig-3 

(0.76-0.510 x iq-3 

(-1.7 + 0.290 x IQ' 3 



v 

(e) Clamped panel; 7 - 2. 

w 

Figure 5.- Continued. 
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Point A 

Point 6 

Point C 

Point D 

Point E 

0.240 

0.1610 

0.0861 

0.0605 

0.0319 

0 

0 

0.010 

0.010 

0 

2.33 

7.0 

8.0 

8.0 

1.0 

2.295 

2.020 

1.691 

1.472 

0.9881 

1.015 

3.46 

4.72 

5.43 

1.012 

-0.67 - 0.41i 




-0.89 - 0.23i 

1.0 




1.0 

-0.61 + 0.321 




-0.48 + 0.085i 

0.20 - 0.211 




0.15 - 0.00351 

-0.061 + 0.082i 

-0.012 + 0.0063i 

(-0.29 + 0.781) x 10‘ 2 

(0.11 - 0.048i) x io-2 

-0.047 + O.Olli 

0.020 - 0.0321 

0.015 - 0.012i 

(0.31 - 0.17i) x 10-2 

(-0.24 + 0.581) x 10' 2 

0.015 - 0.0038i 

(-0.63 + 1.3i) x 10' 2 

-0.035 + 0.022i 

(-0.66 + 0.28i) x 10-2 

(0.24 - 0.076i) x io-2 

(-0.59 + 0.14i) x 10-2 

(0.18 - 0.59i) x 10-2 

0.053 - 0.045i 

(0.81 - 0.35i) x 10-2 

(-0.52 + 0.21i) x 10-2 

(0.21 + 0.066i) x 10-2 

(0.29 + 2.7i) x 10-3 

-0.28 + 0.16i 

-0.019 + 0.0098i 

(0.55 - 0.12i) x 10'2 

(-0.85 + 0.30i) x 10-3 

(0.20 - 1.5i) x 10'3 

1.0 

0.024 - 0.009r 

-0.015 + 0.0085i 

(0.27 - 0.17i) x io-3 

(0.30 + 0.81i) x 10-3 

-0.25 + O.lOi 

-0.17 + 0.090i 

0.013 + 0.00 17 i 

(-0.12 - 0.093i) x io-3 

(-0.30 - 0.531) x 10' 3 

0.031 - 0.0411 

1.0 

-0.15 + 0.052i 

(0.17 - 0.65i) x 10-4 


-0.023 + 0.015i 

-0.16 + 0.057i 

1.0 



(0.61 - 0.98i) x 10-2 

0.013 - 0.0161 

-0.10 + 0.034i 



(-0.64 + 0.48i) x 10’ 2 

-0.014 + 0.0064i 

(0.9 - 0.851) x 10 - 2 



(0.22 - 0.38i) x 10-2 

(0.30 - 0.39i) x iq-2 

(-1.1+0.380x10 
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(f) Clamped panel; ^ = 4. 


Figures.- Concluded. 
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mode- 10 boundary. For g ^ = 0.01, the flutter boundary is formed by portions of the 
mode-2, mode- 12, and mode- 13 stability boundaries. For g i greater than about 0.02, 
the mode- 2 boundary is the flutter boundary in the range of 1/p. shown. For points A 
and E on the highly coupled boundary, the proportion of natural mode 2 is the largest 
present as indicated in the table of figure 5(f). An examination of the generalized coor- 
dinates (the eigenvectors) indicates that the region of the panel undergoing the greatest 
amplitude is located well back in the rear half. For each of the points B, C, and D, one 
mode is predominant in the motion. When these three points were computed, modes 5 to 
16 were used in the analysis. The smallness of the coordinates q^ and q^g are con- 
sidered as an indication that a substantially converged result has been reached. For 
aluminum panels at sea level, the flutter boundary prescribes t = 0.00180 for g. = 0 

f L 1 

and j = 0.00178 for gj = 0.01. The parabolic boundary from the closed-form solu- 
tion falls near the mode-2 boundary and on the conservative side of it. 

It should be pointed out that in the section "Results and Discussion" of reference 1, 
a certain statement lacked a needed qualification. The statement applied to aluminum 
panels at sea level and was "For all the length-width ratios at least up through 4 on 
Fig. 3 (of ref. lj , the flutter frequency was near and usually slightly above the first- 
natural-mode frequency." The statement should have been qualified that for certain 
of the l / w values, the damping of the panel must be equivalent to a value of gj of 
about 0.03 or greater for this generality to be true. 

In the section "Analysis" it was stated that the results presented were obtained on 
the basis of assumed beam modes for both the stream and cross- stream directions. For 
the panel with all edges clamped, the effects of stiffness cross-coupling between modes 
are neglected. This cross- coupling occurs because for clamped edges, the integral over 
the panel of the product of lq times the second term on the left-hand side of equation (3), 

— - — -, is nonzero, where the number of half-waves p in the stream direction for mode 
3x^9y^ 

i plus the corresponding p for mode j is an even number, and the number of half- 
waves q in the cross- stream direction for mode i plus the corresponding q for 
mode j is also an even number (in modes i and j, q = 1 for the present results). 

(For sine-wave mode shapes of simply supported panels, the stiffness cross-coupling 
integral is zero; only the direct term is nonzero where p for mode j equals p for 
mode i.) 

A recent report (ref. 9) gives some quantitative results obtained by including the 
stiffness coupling terms in a Galerkin modal analysis and using aerodynamic forces from 
static strip theory. The dynamic -pressure parameter A increased by 29 percent for 
~ = 4, and by 12 percent for ^ = 2; these increases corresponded to decreases in required 
panel thickness of about 9 and 4 percent, respectively. Unpublished results obtained at 
the Langley Research Center agree with the results cited from reference 9, and also show 
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that with enough modes, the Galerkin solution converges closely to the closed-form result 
from the Kantorovich method of solving the partial differential equation essentially as in 
section 8 C of reference 6, except for the minor difference of using the clamped- end-beam 
fundamental mode for the cross-stream variation instead of the 1.0-minus-cosine varia- 
tion indicated in equation (8.37) of reference 6. 

Thus, for sufficiently high values of 1/ w, the effects of stiffness cross-coupling 
become significant and should be taken into account where the more precise results are 
needed. 


Effect of Mach Number 

Flutter boundaries have been computed for an aluminum panel at sea level with 
— = 2, g^ = 0, all edges clamped, and for several Mach numbers in the range from 1.02 
to 2.0. The results are given in table in in terms of t/l, (/3E/q)^ 3 t, /l, and 
(ME/q 

The thickness ratio t/l required to prevent flutter is virtually constant for Mach 
1.2 to 1.5, is 5 percent higher at Mach 2.0, 6 percent lower at Mach 1.1, and 16 percent 
lower at Mach 1.02. This trend for ^ = 2 in the low supersonic Mach number range is 
in sharp contrast to the result for two-dimensional panels = 0), for which a great 
increase in thickness ratio is predicted for Mach numbers less than about y2. The ratios 
of flutter frequency to first natural frequency fell between 1.0 and 1.1 for these solutions 
for values of g. near zero. Realizable amounts of structural damping have very little 
effect on solution values in these cases. 

The parameter (j3E/q)*/ 3 t/Z is tabulated as a matter of interest because it is in 
rather widespread use, having arisen on the basis of the static aerodynamic approximation. 


TABLE III.- CLAMPED- EDGE ALUMINUM PANELS WITH 
^ = 2 AND g i = 0 AT SEA LEVEL 


M 

t/l 

(/3E/q) 1/3 t /l 

(ME/q) 1/3 t/i 

1.02 

0.00259 

0.150 

0.258 

1.05 

.00259 

.166 

.246 

1.1 

.00292 

.212 

.284 

1.2 

.00311 

.241 

.294 

1.3 

.00312 

.251 

.291 

1.4 

.00311 

.248 

.280 

1.5 

.00311 

.247 

.272 

2.0 

.00327 

.248 

.260 
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For this panel the parameter is nearly constant for M = 1.2 to 2.0, but drops toward 
zero in the manner of as M approaches 1.0. The parameter (ME/q)^^t /l is 

also tabulated as a matter of interest since some investigators have applied piston-theory 
aerodynamics to the analysis of panel flutter. It has the feature of remaining finite as 
M approaches 1.0 for this panel. The values of the two parameters, one with /3 and 
the other with M, approach each other as M increases. 

CONCLUDING REMARKS 

A panel flutter analysis procedure has been developed for which the generalized 
aerodynamic forces are computed by considering the panel to be finely divided into a 
large number of boxes and computing the matrix of aerodynamic influence coefficients 
relating all pairs of boxes. This relatively laborious technique can be resorted to when 
the aerodynamic forces from simpler expressions are either not obtainable or are of 
questionable validity. The technique has particular usefulness in applications to finite 
panels at low supersonic Mach numbers. Flutter solutions are obtained from a modal 
type of analysis, and all the modal quantities are input information. Either experimental 
or analytical mode shapes, frequencies, and mass data can be used. Thus, this type of 
[ analysis can be applied to any essentially flat panel, whether unstressed or stressed (as 
by thermal expansion) or whether of isotropic or anisotropic stiffness, and to a small- 
amplitude flutter superimposed on a buckled deflection. 

Results are presented only for flat unstressed isotropic rectangular panels with 
side edges alined to the airstream direction and with no pressure difference that would 
tend to bulge the panel. Design-type plots of flutter boundaries are presented for Mach 
1.3 for length-width ratios ranging from 0 to 10 for simply supported edges and from 0 
to 4 for clamped edges. Additional information is tabulated for a number of points on the 
flutter boundary of each figure, frequency and flutter-mode-shape information being given 
to aid in a fuller understanding of the results. For length-width ratios of 1.0 and greater, 
flutter boundaries obtained from closed-form (nonmodal) solutions based on static strip- 
theory aerodynamics are included for comparison. A comparison of the dynamic- 
pressure parameter and of thickness required to prevent flutter shows surprisingly good 
agreement with the results of the present method for certain ranges of length-width ratio 
and ratio of air mass to panel mass. These ranges include aluminum panels at sea level, 
having length-width ratios of 2, 4, and 10 with simply supported edges, and length-width 
ratios of 2 and 4 with clamped edges. 

The effect of Mach number variation for a clamped- edge aluminum panel with a 
length-width ratio of 2 at sea level was studied. Three different panel flutter parameters 
are given for Mach numbers in the range from 1.02 to 2.0. With increasing Mach number 
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over this range, the thickness required to prevent flutter increases somewhat irregularly 
by about 20 percent. 

For clamped-edge panels, the stiffness coupling effects from the assumed beam 
modes were neglected. A recently published reference reports that inclusion of these 
coupling effects increased the dynamic-pressure parameter by 12 and 29 percent for 
length-width ratios of 2 and 4, respectively, as determined on the basis of static strip- 
theory aerodynamics in a Galerkin analysis using up to 20 beam modes. These and other 
supporting results indicate that for sufficiently high length-width ratios, the stiffness 
coupling between beam modes should be taken into account for flutter analysis of clamped- 
edge panels. 

Appendix A gives the form of the expressions used for calculating the mode- shape 
information used in obtaining the presented flutter boundaries. The form given was 
adopted because it provides full single-precision accuracy of the mode-shape quantities 
without the need for multiprecision arithmetic for the higher modes. Appendix B pro- 
vides conversion formulas for a number of flutter solution parameters in current use. 
Appendix C describes a way to reduce computing machine time for the large matrix 
multiplications required to obtain the velocity potentials. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., June 19, 1966, 

126-14-02-01-23. 
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NATURAL MODE CHARACTERISTICS OF RECTANGULAR PANELS 


The panel flutter computing program described in the text is designed to accept 
input information for the downwash, that is, the streamwise slopes and modal deflections 
obtained from either experiment or analysis. All the panel flutter boundaries given in 
the present report were obtained with analytically determined downwash input. A 
description of the derivation of the natural mode shapes employed is given. 


For the panel natural vibration modes, it is assumed that separation of the time 
variable and the two space variables is valid; that is, 

Hj(x,y, r) = q-j(r) hj(x,y) = q.e lu,T Xj(p,x) Yj(q,y) (Al) 

where the generalized coordinate qj can be complex to account for leading or lagging 
phase relationships in the general flutter motion, and where Xj(p,x) and Yj(q,y) con- 
tain the x- and y-variations, respectively, and p and q are the numbers of half-waves 
of deflection in the x- and y-directions, respectively. 


The governing differential equation for small- deflection purely flexural vibrations 
of a uniform thin isotropic panel with no in- plane loading and unaffected by any sur- 
rounding air is 


D 


l^j 

8 x 4 


\2 9 4 Hi , , \4 8 4 H j 4 9 2 Hj 

I 4 +m A l ~~T =0 
8 y 4 A 




(A2) 


As done by several previous investigators (see, for example, refs. 10 and 11), tne mode 
shapes Xj and Yj of the plate are approximated by the mode shapes for a uniform 
beam vibrating in flexure. In succeeding sections, certain combinations of edge support 
are considered. These combinations are (1) all edges simply supported; (2) all edges 
clamped; (3) leading edge clamped, trailing edge simply supported, and side edges either 
clamped or simply supported. 


Panels Simply Supported at All Edges 

The mode shape functions can be given in a form that satisfies both the boundary 
conditions and the differential equation (A2) and these functions are 


Xj(p,x) = sin KpX 
K p = Ptt 


(p = 1, 2, 3, . . .)> (A3 a) 
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Yj(q,y) = sin Kqy 

Kq = qir 


(q = 1, 2, 3, . . .)> (A3b) 


The natural resonant frequency of mode j as derived by application of Rayleigh’s 
principle of equal maximum kinetic energy and elastic potential energy is 


,, 2 . D 

i 


far 

(fa ) 4 = V + 2 (wr(V)(-V) + 


(A4) 


where the combination of indexes p and q is appropriate to mode j. 


Panels Clamped on All Edges 

For uniform panels the mode shape functions are approximated by uniform-beam 
mode shapes in both x- and y-directions. More than one form for Xj and Yj is pos- 
sible. The following form was adopted because single-precision (nominally, eight- 
significant-figure) arithmetic operations are sufficient for computing single- precision 
results for the high modes as well as for the low modes: 


Xj(p,x) = ^p[(l - a PjX )e Kp " + (l + Op >x )e p " + 20^ sin K p x - 2 cos K p xj 


-K n x 


(A5a) 


Yj(q,y) = ^ 


(1 - °q,y) eKq ^ + (l + y sin Kqy - 2 cos Kqy 


(A 5b) 


where the amplitude factors A x and Ay are chosen as unity so that 


r x i 


2 dx = A x 2 = 1 


£ ■ 


Ay^ = l 


(A6) 


The characteristic values K p and Kq are the roots of the characteristic equations 

(A7a) 
(A7b) 


cos K p cosh K p = 1 
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from which it is found that 

Kp = - ep (| e p| <<c (A8a) 

Kq = - e q (|e q | « l.o) (A8b) 

The small quantities €p and e q rapidly become even smaller as p and q increase, 
and the labor of determining them is eased by using their asymptotic values for the 
higher values of p and q: 

eq ~ 2(-l) < *e"^ q+ ^ 11 
The quantities ap X and y are 

cosh Kp - cos Kp 
a P’ x sinh Kp - sin Kp 

cosh K q - cos K q 
sinh Kq - sin Kq 

For the evaluation of Xj and Yj from equations (A5a) and (A5b), the values of 
1 - Op X and 1 - Oq y to single-precision accuracy are required; for this evaluation 
their asymptotic values for high p and q are helpful 

1-Q p,x~ -e P (Alla) 

1 - Q q,y~' e q (a lib) 

Values of Kp, Kq, 1 - «p ?X j anc * 1 - f° r values of p and q ranging from 1 

to 18 and for opposite- edge pairs clamped are given in the following table: 


(A9a) 

(A9b) 

(A 10a) 
(A 10b) 
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p or q 

K p or Kq 

X >> 

I O l 

o o 

• • 

tH t-H 

1 

4.7300407 

0.17497779 x 10" 1 

2 

7.8532045 

-.77722999 X 10 -3 

3 

10.995608 

.33551562 X 10“ 4 

4 

14.137165 

-.14498945 X 10“ 5 

5 

17.278760 

.62655621 x 10“ 7 

6 

20.420352 

-.27075949 x 10“ 8 

7 

23.561945 

.11700579 x 10~ 9 

8 

26.703538 

-.50562789 x 10 -11 

9 

29.845130 

.21850161 x 10~ 12 

10 

32.986723 

-.94423106 x 10 -14 

11 

36.128316 

.40803924 X 10" 15 

12 

39.269908 

-.17632974 x 10 -16 

13 

42.411501 

.76198991 x 10“ 18 

14 

45.553093 

-.32928570 X 10~ 19 

15 

48.694686 

.14229725 x 10 -20 

16 

51.836279 

-.61492218 X 10" 22 

17 

54.977871 

.26573197 x 10~ 23 

18 

58.119464 

-.11483320 x 10' 24 | 


For an isotropic material, application of the Rayleigh principle leads to the expres- 
sion for the natural frequency of mode j of a plate 



( K p,q) K P + (w) ^ + 2 (w) “p,x K p( 2 “ K P°fr, x) 0 ^, y Kq ( 2 " Kq^y)! 

where the stiffness coupling effects are neglected. The combination of p and q is 
that appropriate to mode j. 

As was pointed out at the beginning of this section, equations (A5a) and (A5b) give 
approximations for the plate mode shape. The boundary conditions are satisfied but the 
differential equation (A2) is satisfied only approximately. The differential equation for a 
plate could be satisfied with greater accuracy by using a series of terms of the beam- 
mode type. What is done here, as has been done by many investigators, is to use only 
the one predominant term in a series for each mode j and to neglect the other terms on 
the grounds that their effect is minor, not of the first order. Recently published findings 
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of reference 9 regarding the effects on some flutter solutions of including the stiffness 
coupling terms are described for clamped- edge panels in the section "Results and 
Discussion." 


Panels Clamped at Leading Edge and Simply Supported at Trailing Edge 


For a panel clamped at its leading edge and simply supported at its trailing edge, 
the function Xj is again approximated by that for a beam with the same boundary con- 
ditions. Thus, Xj has the same form as in equation (A5a) but the characteristic values 
of Kp are obtained from the characteristic equation 

tan K p = tanh K p (A13) 

the roots of which are 



and for the higher values of p the asymptotes are 



^ ep « 1^ (A14) 


(A 15) 


1 - °p,x ~ - 2ep (A16) 

The quantities Kp and 1 - ap ;X are given in the following table for leading edge 
clamped, trailing edge simply supported, and for values of p from 1 to 6: 


p 

K P 

" °p,x 

1 

3.9266023 

-0.77731188 x 10" 3 

2 

7.0685827 

-.14498977 X 10“ 5 

3 

10.210176 

-.27075950 X 10 -8 

4 

13.351769 

-.50562785 X 10 -11 

5 

16.493361 

-.94423105 X 10“ 14 

6 

19.634954 

-.17632974 X 10' 16 


If the amplitude factor A x is chosen as unity, 



Xj 2 dx = A x 2 


= 1 


The natural frequency of the panels depends on the type of support at the side edges. 
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Simply supported side edges .- For simply supported side edges Yj and Kq are 
given by equation (A3b). There is no stiffness coupling and the natural frequencies are 
given by 



(A17) 


Clamped side edges .- For clamped side edges Yj is given by equation (A5b); Kq, 
by equation (A7b); and a , by equation (A 10b). (The latter two quantities are obtained 
from the table following equation (All).) When stiffness coupling is neglected, the natural 
frequencies are given by 



“\ 

2 

+ 2 (w) ^x 1 ^ 1 " Kp^x^y 1 ^ 2 " V^y) 


(A 18) 


Panels With Various Edge Supports, Orthotropicity, and In- Plane Loadings 

A recent significant contribution to the analysis of vibration of rectangular panels 
is made by reference 12 which treats a variety of edge restraints including simply sup- 
ported, clamped, and elastic restraint against rotation. Account is taken of in-plane com- 
pressive and tensile loads in both directions and of orthotropic stiffness of the panel. 
Analysis is made by the method of Kantorovich, which involves an assumption of the 
modal deflection shape in one direction, a consequent reduction of the partial differential 
equation in two independent space variables to an ordinary differential equation in a single 
independent space variable, and its eigensolution. Solutions that satisfy the plate differ- 
ential equation exactly are obtained if the assumed deflection mode in one direction is 
between simply supported edges. Extensive tables are provided in reference 12 for the 
use of the analyst. 
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RELATIONS AMONG SOME PANEL FLUTTER PARAMETERS 
IN THE LITERATURE 


The flutter boundaries of the present report are presented in terms of a mass 
ratio 1 / /jl plotted against a stiffness parameter cu^^/V. Their relation to some other 
parameters in the panel flutter literature is explained briefly here, specifically for 
isotropic panels only. 

A parameter having rather wide use (in refs. 4 and 5, for example) in a variety of 
forms is a panel flutter dynamic-pressure parameter defined by 


_ pV 2 l 3 _ 2qZ 3 24q(l - v 2 )n$ 

/3D 0D /3E W 


(Bl) 


where q is stream dynamic pressure, and E is Young's modulus of elasticity. The 
relationship of X to the parameters of the present report is 


= (k* Y 4 ! Vm 
1,1 P ( u i l / V f 


(B2) 


and 




is given in appendix A for several edge- support conditions. (Assign 

p,q = 1,1.) The dynamic-pressure parameter is seen in a variety of forms; some of 

these replace /3 by M, drop /3 entirely, drop the factor 24, drop the quantity (l - vty, 

introduce the factor and so forth. Other variations are the cube root and the inverse 
1/3 

cube root X ' . The inverse cube root has the merit that the thickness to prevent flut- 
ter, and thus the weight penalty, are contained linearly. The relationships among such 
variations of X are readily apparent. 

In reference 13 the flutter boundaries are on plots for which the relation of the 
coordinates to X and to the present parameters may not be as clear. The vertical 

)l/3 


coordinate is r 


pc. 




which is equal to 


nate is 


report. 


C ooPq 

Ps(i - v2) 

p 

L E J 

13, T, 

c„, and 


1/2 


which equals 
t/l. 


are 


S’ 


1 "1* 
P V 
and m 


m 

V 1 ^) 2 

/ MN/12 


The horizontal coordi- 


The symbols from ref* 


Mn/12 

j^Jt, respectively, in the present 
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MATRIX OPERATIONS FOR COMPUTING THE VELOCITY 
POTENTIAL MATRIX OF EQUATION (27) 


The practicality of this box method is largely dependent upon the cost of the com- 
puting machine time required. This time is significantly influenced by equation (27). 

Not only should the multiplication be carried out completely within the machine core, but 
also the matrices should occupy as small a portion of the core as possible, consistent 
with the desirably low machine time for the overall program. Thus, any repetition of 
the storage of a given number should be avoided. 

If the |a^ (r , s)j> were indeed a single- column matrix, the rectangular matrix of 
complex downwash ratios would have B g B xs rows and nearly 2B S B XS columns, 
nearly half of the elements would be zero, and all the nonzero elements would be repeated 
from at least Bx S to as much as B S B XS times. But by an arrangement of the sub- 
matrices that can be described as "folded," repetition of the submatrices of the downwash 
is avoided. 


As an illustration, for the simple case of B s = 5 and ®xs ~ 3, let the elements of 
3h-j 

be represented by their indexes ii,v, and let zero 


-fUifhi 

dX V 1 


the downwash matrix 
elements be indicated by a blank. Then that "folded" matrix is as follows: 


4,2 

3,2 

2,2 

1,2 0, 

3,2 

2,2 

1,2 

0,2 

2,2 

1,2 

0,2 


1,2 

0,2 



0,2 





4,1 

3,1 

2,1 

1,1 0,1 

3,1 

2,1 

1,1 

0,1 

2,1 

1,1 

0,1 


1,1 

0,1 



0,1 





4,0 

3,0 

2,0 

1,0 0,0 

3,0 

2,0 

1,0 

0,0 

2,0 

1,0 

0,0 


1,0 

0,0 



0,0 





(Cl) 


If the elements of the matrix A^,(r,s) are represented by their indexes r,s that 
matrix is, in effect 
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The resulting folded b* (m ’ n: 
m,n is then 


0,0 

0,1 

0,2 

1,0 

1,1 

1,2 

2,0 

2,1 

2,2 

3,0 

3,1 

3,2 

4,0 

4,1 

4,2 

0,1 

0,0 

o,i 

1,1 

1,0 

1,1 

2,1 

2,0 

2,1 

3,1 

3,0 

3,1 

4,1 

4,0 

4,1 

0,2 

0,1 

0,0 

1,2 

1,1 

1,0 

2,2 

2,1 

2,0 

3,2 

3,1 

3,0 

4,2 

4,1 

4,0_ 


matrix of equation (27), represented by its indexes 


4,2 

4,1 

4,0 

3,2 

3,1 

3,0 

2,2 

2,1 

2,0 

1,2 

1,1 

1,0 

0,2 

o,i 

0,0 


(C2) 


(C3) 


In practice an actual arrangement of the A^(r,s) matrix elements as shown in 
matrix (C2) with repetition of elements (due to right- left symmetry) would be prohibi- 
tively wasteful of memory-storage space. It is a simple matter to arrange the elements 
in a single one-dimensional column array and select appropriate subcolumns as needed. 
With a little more ingenuity the necessary subcolumns can be used without any storage 
duplication of elements such as that shown in the middle column of the A^(r,s) 


phj , 

matrix (C2). For the same reason, for the + i -^r hj matrix (Cl) the storage space 

need be no larger than the top row (b s times B xs complex numbers). Successive rows 
can be formed in the same space by shifting elements and inserting zeros. Thus, the 
|b*(m,n)j matrix (C3) is obtainable in a systematic and economic manner for each mode. 


For subsequent use in equation (29) the *(m,n)j matrix must be "unfolded" into 
a single column, the subcolumns of which are the columns of the folded matrix (C3). 
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